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A VISCOUS/POTENTIAL FLOW INTERACTION ANALYSIS 
METHOD FOR MULTI-ELEMENT INFINITE SWEPT WINGS 

By 

F. A, Dvorak and F. A. Woodward 
Flow Research, Inc, 


SUMMARY 


An analysis method and computer program have been developed for the 
calculation of the viscosity dependent aerodynamic characteristics of 
multi-element infinite swept wings in incompressible flow. 

The wing configuration consisting at most of a slat, a main element 
and double slotted flap is represented in the method by a large number of 
panels. The inviscid pressure distribution about a given configuration in 
the normal chord direction is determined using a two dimensional potential 
flow program employing a vortex lattice technique. The boundary layer 
development over each individual element of the high lift configuration 
is determined using either integral or finite difference boundary layer 
techniques. 

Once the boundary layer development is known, a source distribution is 
determined as a function of the calculated boundary layer displacement 
thickness and pressure distributions. This source distribution is included 
in the second calculation of the potential^ flow about the configuration, 
and represents the effect of the boundary layer in the modification of the 
potential flow. Once the solution has converged (usually after 2-5 iterations 
between the potential flow and boundary layer calculations) lift, drag, and 
pitching moments can be determined as functions of Reynolds number. 

The new method has a number of features and capabilities which make it 
a unique method at this time. Some of these features include: 

-The inclusion of methods capable of calculating the boundary layer 
development over infinite yawed wings. 

-The inclusion of normal pressure gradient and longitudinal curvature 
terms in the finite difference program. This has led to much 
improved predictions of the performance of multi-element airfoils in 
two dimensions as compared to the predictions of other methods, 
especially in the calculation of profile drag. 

-The use of source distributions rather than the displacement thickness 
directly to represent the effect of the boundary layer on the potential 
flow. This approach is much more efficient than the alternate 
procedure since the influence coefficient matrix representing the 
geometry of the configuration need be inverted only once. Computer 
time expenditures are consequently much less with the new method. 
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-In the future, the effects of tangential injection and suction 
boundary layer control on aerodynamic performance may be calculated 
as well as the effects of fully three dimensional flow. 


The computer program is written in Fortran IV for the CDC 6600 and 
7600 family of computers. The program occupies 100,000 (octal) words of 
storage and operates in the overlay mode. The program has been structured 
in such a way that extension or replacement of individual calculation 
procedures is straightforward. 
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INTRODUCTION 


Background 


The multi-element wing is an essential component of the high-lift systems 
of existing commercial and military aircraft. Historically the design of these 
systems has been dependent upon experimental verification of predicted aero- 
dynamic performance. This approach has been and continues to be a very costly 
and time consuming venture. The advent of high speed computers and of advanced 
numerical methods is however gradually reducing the reliance on the experimental 
method. Calculation methods now exist which permit the solution of many prac- 
tical engineering problems. The prediction of the aerodynamic characteristics 
of two-dimensional lifting multi-element airfoils including the effect of 
viscosity is an important example of this capability. 

The availability of a three-dimensional version of such a method would be 
of considerable value to the designers of modern aircraft high-lift systems, 
particularly with respect to STOL aircraft, where the design problems appear 
to be the most formidable. Trade-off studies could be made for the design and 
analysis of individual components such as the leading edge devices or the 
slotted flaps. A multi-element wing analysis would have other important 
applications , and it is because of the usefulness of such a method that the pro- 
cedure described in this report was developed. 

The method is currently valid for the infinite yawed wing case, but has been 
structured in such a way that at a later date, it can be extended to the fully 
three dimensional case. The addition of viscous effects is accomplished using 
distributed sources determined from the boundary layer calculations. The need 
to add viscous effects is clearly demonstrated by the results shown in Figure 1.1 
(Ref. 3^.). Obviously, the inviscid solution grossly over estimates the performance 
of the airfoil section. 

It was Prandtl who first suggested adding the boundary layer displacement 
thickness to the original geometry to account for the displacement of the inviscid 
flow streamlines by the boundary layer. This approach has since been used suc- 
cessfully by many researchers. A practical computational difficulty arises with 
this approach however, and that is the need in the potential flow calculation to 
re-invert at each pass a large matrix resulting in large computer time expenditures. 
In order to obtain smooth pressure distributions it is also usually necessary to 
smooth the new geometry before each potential flow calculation resulting in 
further computer time expenditures. An alternative procedure stemming from an idea 
first suggested by Preston, Ref. has been successfully adopted in the computer 
program described in this report. Briefly, a source or sink (negative source) 
distribution is determined as a function of the known displacement thickness, 
entrainment rate, and velocity distributions (q = d/ds (P e U e 6*). With the intro- 
duction of a source distribution a new vortex distribution is determined given 


3 



ft 



FIG 1 1 COMPARISON BETWEEN THEORETICAL AND 
EXPERIMENTAL PRESSURE DISTRIBUTION 
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the original geometry, consequently, there is no need to invert the matrix a 
second time. If several iterations between inviscid and viscous flow are 
required, the potential computer time saving is substantial. 


Problem Definition 

The calculation of the potential flow about a multi-element configuration 
represents the first task of any analysis method. Because the analysis is 
limited to infinite swept wings a two-dimensional potential flow method is 
adequate. Future expansion to the fully three-dimensional case suggested, 
however, that any program be written in modular form in order that the two- 
dimensional method could be readily replaced by a three-dimensional method. 

In the two-dimensional case, a mathematical model is required for the flow field 
about a series of arbitrarily shaped bodies in an incompressible, inviscid flow. 

The model must as a further requirement be able to predict the pressure distri- 
bution at selected off-body points above the flap segments. This is necessary 
because downstream of the wing trailing edge the static pressure normal to the 
flap surface is greatly influenced both by the proximity of the flap to the 
wing trailing edge and by the large surface curvature in the flap leading edge 
region. This static pressure variation (see Figure 1,2) has a considerable 
influence on the development of the combined wing wake-flap boundary layer 
downstream of the wing trailing edge. 

With the potential flow field specified, it is necessary to predict the 
boundary layer development over the multi-element configuration. Calculations 
must include stagnation line initial conditions, laminar, transition and turbu- 
lent boundary layer developments and laminar or turbulent separation predictions 
for each element of the infinite swept wing high lift system. The calculations must 
include accurate predictions of boundary layer development in the regions where 
wing or first flap upper surface and cove boundary layers merge with the down- 
stream flap upper surface boundary layer. This requirement is an absolute 
necessity if accurate drag predictions are to be made. Both longitudinal curva- 
ture and normal pressure gradient terms must be included in the governing 
boundary layer equations as each effect has a significant influence on the 
boundary layer development and subsequently on the section drag coefficient. These 
effects are particularly important in the wing trailing edge-flap leading edge 
region. Once the boundary layer development is known, its effect on the external 
flow must be determined. 

A complete analysis program for the aerodynamic characteristics of multi- 
element infinite swept wings is developed by combining the separate potential 
flow and boundary layer calculation procedures. Iteration between the separate 
procedures results in the prediction of viscosity dependent aerodynamic forces. 

The different parts of the flow about a multi-element infinite swept 
wing high lift system consisting of a leading edge slat, the main wing and 
double slotted flaps are shown in Figure 1.3. The different calculation 
schemes that form the elements or modules of the integrated computer program 
are presented in the following sections. 
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FIG. 1.3 FLOW ABOUT A MULTI ELEMENT INFINITE SWEPT WING 
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POTENTIAL FLOW METHOD 


Configuration Definition 


The multi-element airfoil configuration is represented by pairs of 
surface coordinate points. Each element with the exception of the main 
wing, may be specified in its own or a reference coordinate system. The 
main element must be given in the reference coordinate system. Individual 
coordinate systems are related to the reference coordinate system by pivot 
points. The pivot points are prescribed in both the element and reference 
coordinate systems. In order to loft the configuration, element rotation 
angles must also be prescribed. Given the pivot points and rotation angles 
any element may be translated and rotated to the desired location relative 
to the main element . 

Pivot points may be determined based on such requirements as a 
specified slot gap and wing-flap overlap. Leading edge coordinates usually 
provide a convenient element pivot point although in some cases the hinge 
point of a flap on its mechanical track or linkage mechanism gives a ready 
reference point. Figure 2.1 shows a four element configuration in both 
input and lofted positions. 

If the configuration is made up of a main element and one or more 
slotted flaps, then additional analysis is required to determine flap upper 
surface longitudinal radius of curvature for later use in the finite 
difference boundary layer calculation methods. Accurate calculations of 
curvature require the use of very smooth input data. Because of this, it 
was found necessary to use spline functions to represent the surface being 
analyzed. A spline under tension* is first passed through the coordinate 
points representing the flap upper surface. First derivatives dy/dx are 
then determined from the splined curve using analytic expressions. A 
second spline under tension is now used to represent a curve through the 
calculated first derivatives. This spline is likewise differentiated using 
analytic expressions. Once both first and second derivatives of the 
surface are known the radius of curvature can be readily calculated. Figure 
2.2 indicates the success of this technique in relation to known values of 
radius of curvature for the NACA 4412 airfoil upper surface contour. 


* "Splines Under Tension" - a technique developed by Dr. A. Cline of the 
National Center for Atmospheric Research, Boulder, Colorado, for obtaining 
smooth continuous curves from sets of imput coordinate points. 
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FIG. 2.1 MULTI ELEMENT AIR FOIL LOFTING PROCEDURE 
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CHORD WISE DISTANCE X (%) 


FIG. 2.2 COMPARISON OF EXACT & NUMERICAL CALCULATION 
OF SURFACE RADIUS OF CURVATURE 
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Inviscid Method 


The airfoil and associated flap system in their lofted configuration 
is approximated by a large number of planar segments, or panels, with corner 
points located on the actual airfoil or flap surfaces. The geometry of 
a typical two element system is illustrated below: 



Fig. 2.3 Airfoil Geometry Using Planar Panels 

A triangular distribution of vorticity is located on each adjacent 
pair of panels, as shown above. The vortex distribution is identified by 
the index of the common edge, and is given unit magnitude at that point. 
The normal component of velocity induced by the jth vortex distribution at 
the center of panel i is designated the aerodynamic influence coefficient 
a., and is calculated as follows: 


w. . 

ij 



Fig. 2.4 Aerodynamic Influence Coefficients 
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First, the horizontal and vertical components of velocity u.. and w. . 
are obtained by summing the influences of a linearly varying^ortex ^ 
distribution on panel (j-1) having unit value at the trailing edge, and a 
linearly varying vortex distribution on panel j having unit value at the 
leading edge. Formulas for the u and w components induced by these 
vortex distributions in terms of the primed coordinate system associated 
with the influencing panel are given in Appendix I. 


u . . = u ! 1 . cos6 . - - w! 1 . n sin6 . - + u ! , cos6 . - w! . sin6 . (2.1) 

iJ J-1 i,3-l 3-1 iJ 3 1J 3 

w. . = w! 1 . „ cos 6 . - + w! ? . i sin6 . _ + w! . cos6 . + w 1 . . sin6 . (2.2) 

iJ 1,3-1 J-1 i,3-l 3-1 iJ 3 iJ 3 

The normal velocity a^_. is then 


a. . = w. . cos6. - u # sin6. 
iJ iJ i iJ i 


(2.3) 


A series of overlapping triangular vortex distributions are placed on the 
upper and lower surfaces of the airfoil, as indicated: 



Fig. 2.5 Vortex Distribution on Airfoil 


It should be noted that the number of panels on the upper and lower 
surfaces are not necessarily equal. At the leading edge, the vortex 
strengths of the upper and lower vortices are set equal, to insure a smooth 
flow around the leading edge. At the trailing edge the "Kutta" condition 
specifies that the magnitudes of the surface velocity on the upper and lower 
surfaces have a common limit. This implies that the vortex strengths on the 
upper and lower surfaces must be equal and opposite. 
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In the above example, the airfoil has 8 panels on the upper surface and 
7 on the lower, for a total of 15 panels. If the leading edge vortices 1 and 
l 1 are set equal (y^ - y|) and the trailing edge vortices 9 and 16 are set 
equal and opposite ty^ = - y^g) a tota l of 15 unknown vortex strengths 
remain. The unknown vortex strengths are determined by specifying that the 
sum of the induced velocity and the normal component of the free stream velocity 
go to zero at each panel control point. For an airfoil having N panels, the 
total normal velocity at panel i may be written 


N 

n, = sin(a- 6.) + . Z, a . . y . = 0 
i i j=l ij j 


(2.4) 


The first term represents the normal component of a unit free stream 
velocity at the control point of panel i, and the second is the sum of the 
products of the influence coefficients and the N unknown vortex strengths. 
Writing this boundary condition equation for each of the N panels results 
in a linear system of N equations in the N unknown vortex strengths. 

In matrix form. 


- a ll a 12* * * ' * a lN - 
a 21 

a Nl a NN 


V sin <n- fi x ) 


f 1 1 

1 


' r ' 

Y 2 




► { 


, > 

l 

k j 

y n 


sin ( a_(S N ) 


(2.5) 


This matrix equation can then be solved for the vortex strengths. Direct 
inversion is employed for single element airfoils, and either a direct 
method or an iterative procedure (described in Appendix II) can be employed 
for multi-element airfoils. 

Airfoils with blunt trailing edges can be analyzed successfully using 
the Kutta condition that the trailing edge vortex strengths are equal and 
opposite (i.e. y^ = - y^, in Figure 2.5). If the trailing edge closes to 
a point, the strengths of the trailing edge vortices must go to zero, since 
the trailing edge will be a stagnation point in the flow. Although this 
result is given automatically by the solution of the above system of equations, 
it has been found that an alternate formulation of these equations is desirable 
for airfoils with trailing edge closure. In this case the coefficients in 
the last column of the matrix (eqn. 2.5) become very small resulting in a 
poorly conditioned system of equations. 

In the alternate formulation, the Kutta condition is specified by 
setting the strengths of the vortices associated with the trailing ege panels 
equal to zero (i.e., y Q « Y-, - 0 in Figure 2.5). However, this procedure 
eliminates the last column or influence coefficients in Eqn. (2.5) leaving an 
indeterminate system of N equations in N-l unknowns. 
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An additional unknown is provided by adding the influence of a constant 
strength source distribution just inside the airfoil surface. The source 
is distributed on the inner side of each panel used to represent the airfoil. 
The velocities induced by a constant strength source distribution are given 
in Appendix I, and are used to calculate new values for the last column of 
influence coefficients in the boundary condition equations. The unknown 
source strength is added to the remaining N-l unknown vortex strengths to 
give a well conditioned set of equations. It should be noted that the 
unknown source strength is always very close to zero for airfoils with 
trailing edge closure. 


The pressure coefficient at the mid point of panel i is calculated as 
follows: 


where 



W-J 


N 

£ 

3=1 




Y . 
3 


N 

E 

3=1 


w. 


i 3 


Y . 
3 


( 2 . 6 ) 


and u.., w^* are given by Eqns (2.1) and (2.2). The lift and pitching 
moment coefficients are obtained by integrating the pressures around the 
airfoil configuration. 


Viscous/Potential Flow Interaction 

The inviscid flow around an airfoil can be modified to account for viscous 
effects through the addition of the boundary layer displacement thickness 6* 
to the original airfoil geometry. The potential flow method described in the previ- 
ous section can then be used to calculate the flow field about the new geometry. 

A modified pressure field results, which in turn causes a change in the calcu- 
lated boundary layer development. After several iterations both the pressure 
field and boundary layer developments should become convergent. This 
procedure is used in both the Lockheed and McDonnel -Douglas programs (Refs _1^ 
and _3) , and while it would seem at first glance to be a straightforward 
approach, several requirements are necessary to ensure a satisfactory solution. 

These include: 

The necessity to calculate and invert a new influence coefficient 
matrix at each iteration due to the change in resultant geometry. 

The necessity to smooth the geometry each time the displacement 
thickness is added, to ensure a smooth pressure distribution. 
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The necessity during each iteration to modify to a considerable 
extent the calculated pressure and boundary layer developments in the 
trailing edge region in order to ensure a convergent calculation 
procedure. 

An alternate procedure is available which uses the same influence coefficient 
matrix throughout the calculation, and it is this method which is used in the 
present program. 

The effect of the boundary layer on the potential flow is represented by 
a distribution of sources on the panels connecting points of the original air- 
foil surface (Figure 2,3). The strength qj of the source distribution is made 
proportional to the rate of entrainment of mass into the boundary layer (i.e. 
q = (pU e <$*))t Thus the calculated pressure distribution and boundary 

layer displacement thickness developments can be used to determine the source 
strengths for the next iteration of the analysis. The source distribution 
has the effect of modifying the boundary conditions to the original problem 
by altering the right hand side of Eqn 2.5. The influence coefficients a^ 

(or ui j , w-j.) remain unchanged as does the geometry of the configuration 
being analyzed. 

The effect of the source distribution on the boundary conditions is 
determined in the following manner. Consider a panel representing a portion 
of the airfoil geometry; the source strength is known as is the normal velocity 
induced by the source distribution at the boundary point on the panel. This 
normal velocity is the new boundary condition to be satisfied by all sources 
and vortices representing the geometry and the boundary layer effects. However, 

the source distribution of the same panel already satisifies this new boundary 
condition, therefore, the remaining sources and vortices must satisfy the 
boundary condition of tangential flow to the surface. 

Source influence coefficients u sij and w si* are defined as induced 
velocities per unit source strength qj at ^ a corner point where the 

source distribution on a panel is represented by two overlapping triangles. 

This definition is completely analogous to the vortex influence coefficients 

u-l- and w^.. 

J 

The total induced velocities at the i-th boundary point can be described by 


u . 
r 


N 

Z 

j=i 


u. . 

ij 


y j 


N 

Z 

j-1 


J si 


j J 


w . 


N 

= Z 

j-1 


w . 


Y j 


N 

Z w 

j-1 


s ij q j 


(2.7) 


( 2 . 8 ) 


tWith this technique the normal velocity component at the surface n^ and the 
source strength q ' are equal. 
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Because of the introduction of a source distribution, the Kutta condition at 
the trailing edge of the airfoil takes a form different from those used for 
blunt and closed trailing edges in the inviscid flow calculation. With the 
trailing edge being a stagnation point the sum of the vortex and source 
velocities is zero. With the trailing edge sources known the vortex strengths 
Y and can be determined by the condition that the component of velocity 

normal to the trailing edge panel is zero at the trailing edge on the upper 
and lower surfaces. 



Fig. 2.6 Kutta Condition - Modified by Source Distribution 

From the preceeding figure, 

Y u = C—q^ + q u cos6) /sin0 (2.9) 

= (~q u + cos6) / sin6 (2.10) 

It should be noted that the above equations imply equal pressure 
coefficients on the upper and lower surfaces at the trailing edge. Taking 
the difference of the squares of each equation, 

2 


= Cp 1 
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The addition of the source distribution external to the airfoil modifies 
the normal velocity at the control point of panel i. Referring to Eqn. (2.4), 

N N 

q. = n = sin (a- 6. ) + E a., y. + E b.. q. (2.11) 

1 1 1 j=l 3 j=l 13 3 


where 


- w cos6 . - u 


s. . 


sin<5 . 


13 °ij " 'ij 

is the normal velocity induced by the external source on panel j . 

Since the q^. are known, the right hand side of Eqn. (2.5) becomes 

N 


c . =~sin(a-6 .) - E b. ( q . a 

1 1 j=1 13 4 J q i 


( 2 . 12 ) 


Since the vortex strengths at the trailing edge of the upper and lower 
surfaces are also specified in terms of the external source strength by Eqns. 

(2.9 and (2.10), an additional unknown constant source distribution inside 

each airfoil surface is required to solve Eqn (2.5), as described in the previous 

section. 

The advantage in computer time of this procedure results from having to 
calculate the influence coefficients only once. At each successive iteration 
only matrix multiplication is required to determine the new vortex strengths. 

As in the displacement method the effect of the boundary layer corrections tends 
to cause an overshoot or correction in the pressure field solution at each 
iteration. This undesirable feature is avoided and rapid convergence assured 
if the boundary layer development and resultant source distribution is determined 
from a pressure field weighted using fifty percent of the current solution and 
fifty percent of the previous solution. 

Although the precedure does not require extensive smoothing as in the 
displacement methods some limitation on the source strength is required in the 
trailing edge region if rapid convergence is to be achieved. Rapid growth 
of the boundary layer approaching separation in strong adverse pressure gradients 
(typical of configurations at high angles-of-attack) cause abnormally fast growth 
of the displacement thickness, and in turn the source strength. Numerical exper- 
iments indicate that if a limit is placed on the maximum source strength conver- 
gence can occur between two and five iterations. The calculations also indicate 
that this limit is different for slotted airfoil cases, where the boundary layer 
growth on the flap is very much greater than that typical of single element cases. 
More will be said of this limit in a following section. 
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Because the pressure coefficients are determined from the velocities on the 
boundary points of a panel rather than from the vortex strengths at the corner 
points, the trailing edge pressures are not known a priori. Therefore, they are 
calculated by simple linear extrapolation of the pressures from the last two 
boundary points on the upper and lower surfaces respectively. 

A special situation arises if any element of the geometry does not have a 
closed trailing edge. In this case the solution for the inviscid flow about the 
particular element is determined using the Kutta condition y u = - , No 

internal distributed source is required to complete the problem definition. 
Consequently, when boundary layer effects are included in the first iteration 
(with the Kutta condition determined from Eqns. (2.9) and (2.10)), an internal 
source is required to complete the problem definition and it is necessary to 
recalculate the influence coefficient to include the effect of the distributed 
source. Subsequent iterations require only matrix multiplication to obtain 
the vortex strengths. 
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BOUNDARY LAYER CALCULATION METHODS 


The boundary layer development is 'calculated from the stagnation line of 
each element. For the infinite swept wing case two separate calculation pro- 
cedures are used, each for a particular region of the flow. On the upper 
surfaces of the leading edge slat and main wing and for the lower surface of 
every element of the configuration an integral method is used. This method 
is about 100 times faster than the corresponding finite difference method in 
two dimensions. Economy of computer time is essential in an iterative method 
particularly when a multi-element configuration is considered, if the method is 
to be of practical use to the designer. In all cases where passive blowing 
(slotted flaps) or powered blowing is considered, the finite difference method 
is used. This method could be used for the complete boundary layer analysis 
if desired by the user. 

Descriptions of the individual boundary layer and transition analyses are 
presented in the following sections. 


Stagnation Line Initial Conditions 

It has been predicted theoretically by Cumpsty and Head _4 and Bradshaw 3^ 
amongst others that flow along the stangation line of a infinite yawed wing 
appraoches an asymptotic condition. This condition is one where the rate of 
growth of the boundary layer due to frictional forces is balanced by the 
divergence of the flow from the spanwise to the streamwise direction. Cumpsty 
and Head later demonstrated in an experimental study (Ref 6) their earlier 
theoretical prediction. They were able to show that whether the flow is laminar 
or turbulent its integral properties can be determined directly as a function of 
a single non-dimensional parameter C*. The parameter C* (= V 2 /v dU/dx), where 
V is the spanwise velocity, V the kinematic viscosity and dU/dx the chordwise 
velocity gradient at the stagnation line) is a form of Reynolds number which 
correlates well with the streamwise shape factor H, momentum thickness 0 and 
skin friction coefficient Cf/2. The correlations for H and 0 are presented in 
tabular form in Table _1. Initial integral boundary layer parameters are determined 
from the table for the calculated C*. If C* < 1.35 x 10^ the flow is laminar 
otherwise it is turbulent. The appropriate calculation method is then used to 
determine the downstream boundary layer growth (See Figure 3. 1 ) . 

Integral Boundary Layer Methods 


Laminar Method 


A variety of methods exist for the calculation of laminar boundary layer 
developments, the most general of these being based on finite difference methods. 
In the case of the infinite swept wing substantial regions of laminar flow (10% 
chord or more) , are likely only at the lower Reynolds numbers and sweep angles 
and in cases where large flap deflections result in considerable laminar flow 
on the lower surfaces of the flaps. In these instances the effect of the laminar 
flow on the external flow is negligible and the drag contribution very small. 
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Fig. 3.1 Flow Chart for Boundary Layer Calculations 
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Table 1. Stagnation Line Correlation of C* 
with H and Rq.^. 


24 



Because of this, all laminar boundary layer calculations are made using a two- 
dimensional integral approach along external streamlines. The finite difference 
method to be described in a following section can also be used to determine the 
laminar boundary layer development but for the infinite swept wing case as will 
be seen later it appears to be unnecessary for practical calculations. 

The two dimensional method is an adaptation by Curie _7 of a method 
developed by Thwaites _8. In Thwaite's method the momentum integral equation 

de/dx = C f /2 - (H+2)9/U(dU/dx) 3.1 


is written in the form 



d/dx(K/U) = L/U 

3.2 

where 

K = 0 2 /v (dU/dx) 



L = [£-K(H+2) ] 

3.3 


SL = 6/U(8U/3y) y = Q 



Thwaites used exact solutions to a variety of laminar flows to determine 
the relationship between L and K, 

L = 0.45 - 6 K 3.4 

Curie has pointed out that Eqn. 3.4 is not adequate in flows approaching 
separation, and he has suggested an extension or correction to Eqn. 3.4 giving: 

L = 0.45 - 6 K + g(K,y) 3.5 

The parameter y is a function of both the pressure gradient and the curvature 
or second derivative of velocity. 

y = K 2 U (d 2 U/dx 2 )/(dU/dx) 2 3.6 


Curie rewrote Eqn. 3.5 in the form 

L = F q (K) - yG 0 (K) 3.7 

where F 0 and G Q are universal functions determined from a series of exact 
solutions to laminar flows in the same way as did Thwaites for Eqn. 3.4. 

After substitution of Eqn. 3-5 into Eqn. 3.2 and with subsequent integration, 
the result can be rearranged in the form 

e 2 = 0.45V/U 6 J (1 + 2 . 22g) u 5 dx. 3.8 
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This equation is conveniently solved by iteration, g initially equal to 
zero- With values of K and y determined in the first iteration, a second 
iteration is carried out using Eqn. 3.7. At each step in the calculation 
the local skin friction coefficient, C f /2 and the shape factor H can be 
calculated using Eqn 3.3. The local skin friction coefficient has been 
defined as 

C f = (y/pGU ) a 3.9 

where & in Eqn. 3.3 is determined in a similar manner to L from a series of 
known solutions to give 

l 2 = Fjl (K) - y G x (K) 3.10 

The functions Fq, Fj_, Gq and G-^ are tabulated in the computer program. 
Calculations begin at the stagnation line in the swept wing case, with the 
initial momentum thickness 0 given as a function of C*. If the flow is two 
dimensional K takes a initial value K = 0.0855 at the stagnation point from 
which the initial momentum thickness 0 O is 

0 q = (. 0855 v/dU/dx) 1/2 3.11 

The calculation proceeds either to laminar separation or to the end of 
the airfoil whichever occurs first. The calculated boundary layer development 
is then interrogated to determine if transition, laminar separation or forced 
transition (boundary layer tripping) has taken place. If any of these phenomena 
have occurred the downstream flow is assumed to be turbulent. 

Turbulent Method 

Several methods have been developed for the calculation of infinite swept 
wing three dimensional boundary layers. Among the more useful of the methods 
are those by Cumpsty and Head % Nash If), and Bradshaw ^5. Nash’s method (a 
finite difference procedure) is also applicable to fully three dimensional 
boundary layers but in the words of the originator is cumbersome and inflexible 
when applied to complex geometries. In practical cases the methods of Cumpsty 
and Head and of Bradshaw appear to give similar results, with Cumpsty and Heads 
method having a considerable advantage both in speed and convenience. Because 
of this, their method was chosen for use in the viscous/potential flow interaction 
program. 

In developing their method, Cumpsty and Head chose an orthogonal curvilinear 
system of coordinates based on the projections of external streamlines on the 
surface. In this system streamwise turbulent boundary layer velocity profiles 

resemble very closely two-dimensional profiles. When the streamwise profiles are 
known the cross-flow velocity profiles can be calculated as functions of the 
streamwise profiles and the angle between the surface streamline and the projection 
of the external streamline on the surface (angle g) . 
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Cumpsty and Head wrote the streamwise and cross flow equations in integral 
equation form as follows; 

Streamwise Momentum Equation 


30 


11 


30 


3s 


- K 


12 


3s 


1 3U 

FTT (2 + ® 6 n 


3 o i 

15 <e u -9 22 )K - Cf /2 (3.12) 


Cross Flow Momentum Equation 


^21 k !! 22 + J _!1 

3s K 3s U 3s 
s 


„ 3U aa 

5 e 21 - f 17 (6 11 (H+1) +9 22 } - 2K 8T e 21 

s 


(3.13) 


where k= tan a. 

Equations (3.12) and (3.13) contain several unknowns and before turbulent 
boundary layer predictions can be made further relationships are required 
between the streamwise and cross flow momentum thicknesses, the streamwise shape 
factor H and the streamwise skin friction coefficient C^. 

Entrainment Equations 


The fact that the streamwise velocity profiles are similar to two 
dimensional velocity profiles led Cumpsty and Head to assume that the rate of 
entrainment along a streamline in the infinite swept wing case could be deter- 
mined using relationships developed for two-dimensional flow. This is a 
credible assumption since the entrainment process is a function of the velocity 
defect in the outer part of the boundary layer, a region where the streamwise 
and two dimensional profiles are expected to agree most closely. Cumpsty and 
Heads entrainment equation takes the form 


3(6-6*) 36* 

K 3T 


F (H,) + <«-«*)« |f -i 


3U 

s 

3s 


(3.14) 


where = (5-6*) /0 

The function F(H^) is taken in the form presented by Head 11. 


H ! = FO^) 


(3.15) 
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Likewise the expression relating to the more usual shape factor H is also 
given by Head. 

H = G(H) (3-16) 


The functions F and G in Eqns 3.15 and 3.16 can be analytically defined 
as follows: 

F( Hl ) = exp [-3.512 - 0.617 In (Hj-3)]. (3.17) 

G(H) = 3.3 + exp [0.4667 - 2.722 In (H-0.6798) (3.18) 

for H < 1.6 or 

G(H) = 3.3 + exp [0.4383 - 3.064 In (H-0.6798) (3.19) 

for H > 1.6. 


Streamwise Velocity Profiles 


Cumpsty and Head demonstrated that the streamwise turbulent boundary layer 
velocity profiles could be represented quite accurately by the two dimensional 
velocity profile family derived by Thompson 12^ The law of the wall - law of 
the wake velocity profile family of Coles L3 gives results which are in good 
agreement with Thompsons profiles and could easily be used as an alternate 
approach. 


Cross Flow Profiles 


The cross flow profiles have been specified by the simple relationship 
between streamwise and cross flow velocities suggested by Mager 14, 


v/u = (1 - ? /6) 2 


tang 


(3.20) 


where 3 is the angle between the surface streamline (resulting skin friction 
direction) and the projection of the external streamline on the surface, £ is 
the direction normal to the surface, and u and v are the streamwise and cross 
flow velocities. 
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Cross Flow Thicknesses 


The cross flow thicknesses have been defined using a power law velocity 
profile instead of one of the more complicated two parameter profile relation- 
ships, in the streamwise direction. This approach greatly simpliried the 
definition of the cross flow thicknesses without any great loss in accuracy. 
The use or the power law relationship 


u/U 


(5/6) (H 1)/2 


( 3 . 21 ) 


in Equation (3.20) gives the cross flow thicknesses as defined by: 


5* - e u D(H) ta„6 ; D(H) . -(H + ^ ^ 


9 12 = 6 11 J(H) ta " B * - J(H > - - <H + V ll-j5r + iS2 + S3-H*5 


9 2i ' e u E(H) tanB : E<H) - - (H + V Is - sr + i«l 


(3.22) 


0 22 6 11 tanB 5 + IH H+l + H+2 “ H+3 + H+4 


1 4 L 6 4 


Skin Friction Coefficient 


The streamwise skin friction coefficient is determined using Thompson’s 
two parameters skin friction law although here again Cole’s skin friction law 
could also have been used. The relationship is of the torm 


C f = f(H, R 0 ) 


and is given as C f = exp (AH + B) 


(3.23) 


where 


A = .01952 - . 3868Z + .02834Z - .0007Z J 
B = .19151 - . 8349Z + .06259Z 2 - .001953Z' 
Z = In Re u 


The cross flow skin friction coefficient Cf 2 is then determined from C 
Cf = Cf^ tanB and the resultant skin friction coefficient as C f 



/cosB. 
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Calculation Procedure 


With initial values of 0^ an d H known either from the laminar boundary 
layer calculation or the stagnation line initial conditions Equations 3.12, 3.13 
and 3.14 are integrated using standard integration procedures. The parameters 
0, -j , H and 3 in conjunction with the skin friction coefficient, and the cross 
flow thicknesses are determined along streamlines as functions of the known 
pressure distribution. 

The streamwise free stream velocity U and angle a are determined as shown 
in the sketch, also shown is the angle 3, %he angle between the projection of 
the external streamline on the surface, and the resultant skin friction direction. 



The sum of the two angles a and 3 is continuously monitored during the 
calculation. If this sum reaches 90 the flow is completely spanwise and 
by definition turbulent separation has occurred. The calculation is stopped 
at this point . 

Boundary Layer Transition and Laminar Separation 

Boundary layer transition is a very complex phenomenon and to date no 
reliable theoretical method has been developed for its prediction. Reynolds 
number is a controlling parameter, but it has been shown that the Reynolds 
number at transition can be increased a considerable amount by careful 
elimination of disturbances. At very low Reynolds numbers, laminar boundary 
layers are stable to small disturbances, however, at higher Reynolds numbers 
the boundary layer is unstable, and small disturbances can be amplified. 
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Amplification of these disturbances cause the flow to become turbulent. The 
point at which flow break down occurs depends on the strength and dominant 
frequency of the initial disturbance. Disturbances may be due to freestream 
turbulence, surface roughness, noise or vibration of the surface. As there 
is no detailed analysis of the transition process, transition prediction is 
accomplished by means of empirical correlations. Granville 15^ has developed 
a procedure based on the determination of the neutral stablility point and 
the transition point. The neutral stability point is defined as that point 
downstream of which small disturbances are amplified within the boundary 
layer. It is this amplification of small disturbances that ultimately leads 
to transition. The neutral stability point is reached when the Reynolds 
number based on the local momentum thickness and the local flow properties 
attains some critical value, Rq. . : .S chi ich ting and Ulrich (Ref. 16) have 
show^that R0. can be correlate! with the local pressure gradient parameter 
K = 9 /v(dU/dsy. Correlations by Smith _T7 and others have been reduced 
to analytical form as follows: 

Instability Curves 

K = -0.4709 + 0.11066 In R A 

o 6 

-0.0058591 In R 0 (3.24) 

for 0 < R 0ins <650 


and K = 0.69412 - 0.23992 In R ft 

2 6 
+ 0.0205 In R 0 

for 650 < Ra. < 10,000. 
D ins — 


(3.25) 


If for a given R 0 the pressure gradient parameter K as calculated by 
Eqn. 3.24 or 3.25 is greater than that determined by the boundary layer 
development the flow has passed from a stable to an unstable region. Once 
the flow passes into the unstable region, the transition process begins, 
and Granville has been able to show that a correlation similar to the instability 
process can be used to determine the transition point. 

Granville formed an average pressure gradient parameter K defined as 


K 



trans 

K ds 


ins 


s 


trans 


s . 
ins 


( 3 . 26 ) 
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which correlated reasonably well with the momentum thickness Reynolds number 
at transition This correlation is presented in analytical form as 

follows: trans 


Transition Curves 

K = -0.0925 + 7.0 x 10 _5 R e 

f ° r 0 " ^ trans ± 750 ’ 

K = -0.12571 + 1.14286 x l(f 4 Re 

for 750 < Rq_ < 1100 , 
y trans — 


and 


K = 1.59381 - 0.45543 In R Q + 0.032534 In 2 R 0 

for 1100 < Rq < 3000. 

D trans — 


(3.27) 


(3.28) 


(3.29) 


When the K calculated by one of the above expressions for a given Rq is 
greater than the value determined from the boundary layer development, 
transition is predicted. 

With transition predicted, initial values of the momentum thickness 0 
and the shape factor H are required to start the turbulent boundary layer 
calculation. Because the boundary layer growth is continuous the momentum 
thickness at transition is used as the initial turbulent momentum thickness. 
Since the shape factor varies from values greater than 2.0 to less than 1.5 
through the transition region an empirical expression is used to determine the. 
initial turbulent shape factor. The empirical relation between H and Rq 
was determined from data obtained by Coles 18_: traub 


H 


t 


1.4754 

k°®10^ trans 


+ 0.9698 


(3.30) 


In many cases the pressure gradient is of sufficient strength to 
separate the laminar boundary layer prior to transition. Except in extreme 
cases the boundary layer will then reattach, usually as a turbulent boundary 
layer. Only recently have researchers been able to analyze this phenomenon 
(Ref. 19) and as yet the procedure is extremely complicated and cumbersome, 
consequently empirical relationships are still required. From the measurements 
of Gaster (Ref. 20), and others a correlation is formed which is capable of 
predicting both the occurrence of separation and later the reattachment as a 
turbulent boundary layer or the catostrophic separation. 
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The correlation is of the form: 


and 


K = 0.027 - 0.0007575 Rq - 0.000001157 Rq 2 (3.31) 

for Rq > 125 

K = -0.09 

for R 0 < 125. (3.32) 


If K becomes less than -0.09 separation occurs, and if Rq is less 
than 125 the boundary layer is not able to re-attach. However, if Rq is 
greater than 125 the value of K determined by the boundary layer development 
must be less than that calculated by Eqn. 3.31 before separation without 
reattachment is predicted. If reattachment is predicted, the turbulent 
boundary layer calculation is initiated using the momentum thickness calculated 
at the separation point. 
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Finite Difference Boundary Layer Method 

On elements of a high lift configuration where the flow field is 
particularly complex, such as in the region where the wing wake mixes or 
interacts with the flap upper surface, integral boundary layer methods are 
not capable of completely analyzing the flow, A more satisfactory. method can be 
developed using finite difference methods. Such a method is described in the 
following paragraphs. 

Governing Equations 

The governing equations of mean motion for three-dimensional incompressible 
flow in a general system of curvilinear orthogonal coordinates are: 
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(V. T ) 

ph 3 9y p y 


(3.35) 


Continuity Equation 


■k (h 3 pu) + h (h l h 3 pv) + -k (h l pW) 


= 0 


(3.36) 
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The shear stress terms in Eqs. (3.33) and (3.35) are: 
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(3.38) 


If a practical three-dimensional high-lift system is considered where it can 
be assumed that curvature effects in the spanwise direction are negligible 
compared to the normal chord direction, a surface coordinate system can be 
employed, where 


x - a 

h 1 = 1 + ky 

k = f(x) 

y = 6 

h 2 = 1 


z = Y 

h 3 = 1 



where k is the longitudinal surface curvature. 

The equations can then be written in the following form: 


u 8u jhi dxx uvk _ 1 3P _3_ /Bu \ v k Bu 

X * 1+ky Bx V By W Bz 1+ky p(l+ky) Bx V By \3y/ 1+ky By 



u 2 k = 1 3P 
^ ' 1+ky p 3y 


(3.39) 

(3.40) 
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z : 


(3.41) 


u 3w 3w 3w 1 3P . 3 /3w \ vk 3w 

1+ky 3x V 3y W 3z p 3z V 3y \ 3y / 1+ky 3y 


Continuity 


(pu) + [ (l+ky)pv] + [ (l+ky)pw] = 0 (3.42) 


Equations (3.39) to (3.42) represent the laminar boundary layer equations in 
incompressible flow. The corresponding equations in turbulent flow are: 
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(3.45) 


The continuity equation is unchanged. The terms u'v* and v T w T 
represent the Reynolds stre sses in the normal chord and spanwise directions. 

The shear stress term u T v T is represented by the expression 

- — i — p = v (3u/3y - uk/l+ky) (3.46) 

u v t x 

where v t is the eddy viscosity, which in this case is determined using a two 
dimensional model. Then, if the shear stress vector is considered to be aligned 
with the rate of strain vector (Nash and Patel (21)), the eddy viscosity in the 
" z " direction may be determined from the expression 
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or the equivalent form 


= v t (3.48) 

x z 


Eddy Viscosity Model 

The eddy viscosity model used in the calculation procedure is a modification 
of one developed for two-dimensional turbulent boundary layers and wall jets 
over curved surfaces (Ref. 22). The basic two layer model consists of inner and 
outer regions. The inner region profile is calculated using the modified 
Van Driest relation 


where 


v t (y) = l 2 


1 - exp 
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3u 

dy 


(3.49) 


t = 0.435y 

1/2 

u x = (T w /p) 

A = f (dP/dx, vw) 


The outer region eddy viscosity profile is determined from the Eddy 
Reynolds number (uj o/v t ) and the intermittency function (y(y)). formulations 
described in Reference 22 . These functions are combined to give: 


v t (y) = (ujO/C) y(y) 


(3.50) 


where the velocity scale and the length scale a (standard deviation of 

the intermittency function y) are known functions of the shape factor H and 
the displacement thickness 6* for conventional boundary layers . 
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The extension of the eddy viscosity model to the mixed flow case 
associated with slotted flaps is based on observations made of the develop- 
ment of walls jets exhibiting both velocity maxima and minima in the profile 
(see sketch) . 


In the region below the velocity maxima 
(Region A) the conventional eddy viscosity 
model is used to describe the flow. The 
outer region of the profile (Region C) 
represents the remnant of the upstream 
boundary layer having a large value of 
the shape factor H. Measurements of the 
standard deviation of the intermittency o 
indicate that an asymptotic value is 
approached at high values of H. Similar 
behavior is observed for the velocity 
defect Ud. By employing the asymptotic 
values of er and in the Eddy Reynolds 

number relation u d a / v t» a layer is 

established which when joined to the conven- 
tional two layer eddy viscosity profile 
provides a completed eddy viscosity model. 


Wall jets exhibiting only a maximum 
in velocity have been studied by many 
researchers. Measurements of the stand- 
ard deviation a of the intermittency 
function y for these flows allow 
predictions to be made of the eddy 
viscosity in the outer region of the 
profile (Region B of sketch) . The 
velocity defect in this case is 

simply U max - U e . Region A is again 
described using the conventional boundary 
layer model for the eddy viscosity. 



38 



The slotted airfoil (passive blowing) case results in velocity profiles 
which are very similar to those for tangential injection; profiles with 
velocity maxima in the case of the leading edge slotted slat, and profiles with 
velocity maxima and minima for slotted flaps. As a consequence it was assumed 
that the same eddy viscosity approach as developed for tangential injection 
could be applied to the slotted airfoil case. The initial velocity profile for 
the slotted case is shown in Figure 3.2. Two factors make the problem slightly 
different from the tangential injection case. These are: 

(i) With slotted configurations, the persistence of the potential 
core , and 

(ii) on the flap surface the possibility of considerable laminar flow 
at least to the suction peak. 

Consequently the flow may consist of a laminar boundary layer, above which is 
a potential core. Above the potential core is the remnant of the cove and 
upper surface turbulent boundary layers of the wing or preceding flap segment. 

To account for the presence of the laminar boundary layer and the potential 
core, the eddy viscosity is set to zero in these regions. In the potential 
core region the remaining viscous terms are negligible in comparison with the 
inertia terms, resulting in a form of Bernoullis 1 equation. Three different 
eddy viscosity distributions are possible depending on the flow regime 
(Figure 3.3). 

The inclusion of curvature terms is essential to the success of the cal- 
culation method as is the necessity of including the static pressure variation 
in the direction normal to the airfoil surface. In regions away from the wing 
trailing edge - flap leading edge Eq. (3.44) is adequate; however, if the afore- 
mentioned region is of interest, then a pressure field P(x,y) in the two- 
dimensional and infinite swept wing case or P(x,y,z) in the full three-dimensional 
case must be prescribed. The pressure field above the individual flap surfaces 
is determined directly from the known induced velocity field. 

Once the eddy viscosity distribution, the surface curvature and the pressure 
field P(x,y) are known, Eqns. 3.43, 3.44 and 3.45 can be solved in conjunction 
with the continuity equation. In the present calculation all spanwise gradients 
have been neglected. The resulting equations are solved using a modification of 
the Crank - Nicholson procedure (23) first described in Ref. 24. 


The initial velocity profile is determined by combining known or assumed 
velocity distributions in the wing trailing edge - potential core region. 

The integral method C.IBL) is used to calculate the boundary layer development 
to transition or to some point on the flap surface downstream of the wing trailing 
edge if transition takes place in the slot region. The calculated integral 
parameters at the slot exit are then used to determine the laminar or turbulent 
boundary layer velocity distribution and thickness. If the flow is laminar, the 
laminar boundary layer profile on the flap upper surface is represented by a 
Pohlhausen polynomial. If the flow is turbulent, Thompson’s velocity profile 
family is used to calculate the velocity distribution. 
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FIG. 3.2 INITIAL VELOCITY DISTRIBUTION FOR A SLOTTED 
AIRFOIL CONFIGURATION 
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FIG. 3.3 POSSIBLE EDDY VISCOSITY PROFILES ON A 
SLOTTED FLAP 
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The potential core velocity distribution is determined from the calculated 

off-body pressures (P(x,y)j, The wing lower surface boundary layer profile 

is represented by a power law profile, while the wing upper surface trailing 

edge velocity profile is determined using Thompson *s velocity profile family 

(12) and the values of H, R and at the wing trailing edge. 

0r 


Aerodynamic Forces 

The aerodynamic lift coefficient for a given configuration can be 
determined in several ways. For closed trailing edge airfoils, the most 
accurate procedure involves summing the individual vortex sheet strengths. 


from which 



2r 
u c 


oo 


(3.51) 


(3.52) 


where 


T * circulation about airfoil 
y » vortex strength of i^ singularity 

U - freestream velocity 

00 

C = reference chord. 


When the airfoil trailing edge remains open Eqn. 3.51 does not necessarily 
give the correct circulation even though the vortex distribution is a valid 
solution for the given boundary conditions. The pressure distribution deter- 
mined from the vortex distribution is in very poor agreement with experiment 
(see Figure 3.4). Consequently it has been found that more satisfactory 
pressure distributions can be determined from the expression 


C 

Pi 



(3.53) 


The lift is then determined by integration of the pressure coefficients 
about the airfoil. 
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The flow model used to represent open trailing edge airfoils is being 
reviewed. It is possible that with the Kutta condition y = - Yi a small 
amount of circulation inside the airfoil gives rise to non-zero tangential 
velocities inside the airfoil and consequently vorticity strength which is 
greater than if the flow inside the airfoil was stagnant. 

Pitching moment characteristics are determined as follows: It is 

presumed that incremental resultant pressure forces act at the center of 
pressure of a small panel of a prescribed length. Consequently, the pressure 
force times the moment arm to some reference point gives the increment in 
pitching moment for that point. The sum of the incremental pitching moments 
for each calculated pressure gives the pitching moment for the system. 

The profile drag is determined for a streamwise section of the infinite * 
span wing by using the Squire and Young drag formula (25) . 



< U T.E. /U »> 


H T.E.+5 

2 


(3.54) 


The streamwise momentum thickness 0 , velocity U and shape factor H 
are used in Eqn. 3.54. The skin friction S drag coefficient is determined b 
by the summation of the local skin friction forces in the drag direction. 
Pressure drag is determined by taking the difference between profile drag 
and friction drag. 


CALCULATION PROCEDURES 

All of the calculation methods, potential flow and boundary layer, are 
incorporated into a single computer program. The calculation sequence is 
outlined below: 

(i) The potential flow pressure field is computed for a multi- 
element infinite swept wing configuration (consisting of up to four elements, 
a leading edge slat, the main airfoil, and a double-slotted flap). 

(ii) The boundary layer properties are then computed for each element 
of the configuration as a function of the potential flow pressure distribu- 
tion. Included in these calculations are the locations of transition or 
laminar separation and turbulent separation, if present. 

(iii) Source distributions are determined to represent the displacement 
effects of the boundary layer on each element and of the wing wake-flap 
boundary layer interaction. 
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(iv) A new potential flow solution is then computed taking into account 
the source distribution computed in step (iii) above. 

Steps (ii) through (iv) are^ repeated until convergence (based on the pressure 
distribution and lift coefficient) is achieved, or until the case is abandoned 
for reasons such as large separation zones. Lift, drag and pitching moment 
coefficients are then calculated for the given configuration. The approach 
is illustrated in Figure 4.1. 

The actual program overlay structure is given in Figure 4.2. The main 
supervisor program has been called VIP (for viscous/potential flow interaction). 
This program directs the overall flow of the calculation. The other programs 
include POTFLOW (potential flow), IBL (integral boundary layer method, 

INSPAN (infinite span finite difference method, and FELDPT (field point 
calculation for off-body pressures) . 


CALCULATIONS AND DISCUSSION OF RESULTS 

The ultimate test of any analysis method is in how well does it predict 
actual aerodynamic performance. This can be determined in the case of potential 
flow methods by comparison with exact solutions; for boundary layer methods, 
the usual recourse however, is comparison with experiment. Evaluation of the 
overall viscous/potential flow interaction analysis can also be made only 
through comparison with experiment. The comparisons that follow represent a 
cross-section of the possible configurations that can be treated by the analysis 
method . 

The potential flow method developed as part of the contract effort has 
been compared with several exact' potential flow analyses. Of considerable 
interest is the comparison for the highly cambered Karman-Tref f tz airfoil 
shown in Figure 5.1. Hess (27) has used this case to demonstrate the degree 
of agreement between his new method and other classical methods. It is there- 
fore encouraging to note that our analysis is in almost total agreement with the 
exact case in comparison with other methods. A second comparison with an 
exact solution is for the two element slotted flap airfoil configuration of 
Williams (28). Here again agreement between the numerical approach and the 
exact solution is excellent (Figure 5.2). 

Two calculations have been included to demonstrate some of the capability 
of the finite difference boundary layer method in two -dimens ions (see Ref. 22). 
The effect of longitudinal surface curvature on the boundary layer development 
is shown in Figure 5.3. It will be seen that when curvature effects are ignored 
the calculation is in poor agreement with the data. The calculations shown in 
Figure 5.4 demonstrate that velocity profiles typical of those found on the 
upper surfaces of slotted or blown flaps can be predicted quite accurately. 
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Overlay (0,0) 



FIG. 4.2 VISCOUS/POTENTIAL FLOW PROGRAM OVERLAY STRUCTURE 
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FIG. 5.2 COMPARISON BETWEEN NUMERICAL AND 
EXACT POTENTIAL FLOW SOLUTIONS 
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FIG. 5.3 COMPARISON OF CALCULATED AND 

MEASURED TURBULENT BOUNDARY LAYER 



0 A □ O EXPERIMENTAL DATA 



FIG. 5.4 COMPARISON OF CALCULATED AND MEASURED 
VELOCITY PROFILE DEVELOPMENTS 
DOWNSTREAM OF A BLOWING SLOT 
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Calculations were made for a series of angles-of-attack for the new 
NASA GA(w)-l airfoil. Comparisons were made with the NASA measurements and 
with calculations made by Morgan (29) using the Lockheed program. The 
results are shown on Figures 5.5, 5.6 and 5.7. Both methods give excellent 
agreement with experiment in comparison with measured lift and pitching moment 
coefficients (Figure 5.5). At the higher angles-of-attack the source method 
appears to be in better agreement with experiment. In all cases however the 
Lockheed program is in slightly better agreement with experiments in the 
trailing edge region in the prediction of pressure coefficients. In the 
present program (VIP) the pressure coefficients are calculated on the original 
airfoil surface, and it is believed that if the pressures are determined 
at off-body points defined by the displacement thickness that improved agreement 
with experimental pressures will result. This procedure will be tried at a 
later date. Measured and calculated drag polars are shown in Figure 5.6. 

No allowance has been made for trip drag or separation effects in the calcu- 
lations although these are present in the measurements. Calculated and measured 
pressure distributions are compared in Figure 5.7. Currently neither theo- 
retical approach is capable of predicting the effects of separation present in 
the measurements. A further set of calculations were made for the NACA 23012 
airfoil. Comparisons between theory and experiment for lift coefficient versus 
angle-of-attack (Figure 5.8) and lift versus drag (Figure 5.9) are in good 
agreement . 

The multi-element prediction capability of the program is demonstrated in 
Figures 5.10 through 5.14. The first case considered is that of the NACA 
23012 airfoil with a 25 percent chord slotted flap (Ref. 30). As shown in 
Figure 5.10 the present method is in better agreement with experiment than 
is the Lockheed program. It is believed that the use of a finite difference 
boundary layer method including the effects of curvature and normal pressure 
gradients results in an improved physical representation of the flow in the 
wing trailing edge-flap upper surface region. This in turn results in an 
improved prediction of the circulation about the complete configuration when 
viscous effects are included. Similar conclusions can be drawn from the results 
of Figure 5.11 for the NACA 23012 configuration having a leading edge slot 
and a slotted flap (Ref. 31). 

The NACA 64AQ10 airfoil with, leading edge slot and a double slotted flap 
(Ref. 32) is considered in Figure 5.12, Also shown for comparison are the 
results for the same configuration from the Lockheed program. The comparison 
is similar to that of Figure 5.11 in that the greatest difference between 
the two results is in the prediction of the pressure distribution for the 
main element. It is believed that the difference is a result of the way in 
which the two programs treat the mixing between the main element and the 
double slotted flaps. 
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a , deg. 

FIG. 5.5 LIFT AND MOMENT COEFFICIENTS FOR 
NASA GA (W) -1 AIRFOIL 







FIG. 5.7 PRESSURE DISTRIBUTION NASA GA (W) -1 
AIRFOIL (12.04° ) 
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FIG. 5.9 DRAG POLAR FOR NACA 23012 AIRFOIL 
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O NACA 23012 WITH L.E. SLAT 


AND SLOTTED FLAP CL = 2.08 



FIG. 5.1 1 COMPARISON OF MEASURED AND PREDICTED 
PRESSURE DISTRIBUTIONS FOR NACA 23012 
AIRFOIL WITH L.E. SLAT AND 25% SLOTTED FLAP 






A second and considerably more challenging case has been investigated 
using the same basic four element geometry, In this example the slot and 
double^slotted flaps are positioned in a landing configuration ( ^ s ^ at = 

- 26,1 , <5^ _« 52,7°). Two problems have been encountered at s a 

this time an§^they will require further study. The first problem centers on 
the iterative procedure for solution of the influence coefficient matrix. 

Each element is initially analyzed in isolation with respect to its neigh- 
boring components. Interference effects are then determined and the analysis 
continued until convergence is achieved. The problem arises when interference 
effects are added to the slot loading (Initially strongly negatively loaded 
because of the negative angle-of-attack) , and change the loading to strongly 
positive lift. The result is a divergent solution. Another possible difficulty 
with the configuration analysis may be with the method of solution. Linearly 
varying vorticity methods do not have a strong diagonally dependent matrix as 
in the case of constant source panel methods. Iterative methods of solution 
rely to a considerable extent on the strong diagonal for their success. More 
work is definitely needed in the area of fast reliable solution techniques. 

When the problem with the iterative solution was encountered, the direct 
technique was used to obtain the inviscid pressure distribution. As a result 
of a very high suction peak in the trailing edge region of the upper surface 
of the main element (probably due to the very high camber effect resulting 
from the highly deflected flaps) the starting velocity profile to the finite 
difference program had unacceptably high velocities. More work is needed to 
resolve this problem. 

The RAE 2815 configuration tested by Foster (33) consisting of a main 
element and a single slotted flap has been considered in Figures 5.13 and 5.14. 
The comparisons in Figure 5.13 include measurements and calculations for two 
configurations. The spike in pressure measured by Foster is not duplicated 
by NASA Ames, nor is it reproduced in the calculated pressure distributions. 

This has a marked effect on the velocity profiles shown in Figure 5.14. 

Although the initial velocity profiles are in reasonable agreement, the 
different pressure gradient conditions experienced by the measured and 
calculated boundary layer developments results in quite different downstream 
profiles. It is interesting to note, however, that the aerodynamic load 
comparisons shown in Table 2 are in generally good agreement. Comparisons are 
also made for the RAE 2815 configuration having a drooped leading edge, and a 
slotted flap deflected 30 degrees. All comparisons were at 9 angle-of-attack 
and at a Reynolds number of 3.8 million. 
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O FOSTER 


2.5% GAP, 3.1% O.L. 


□ NASA AMES 1.5% GAP, 1.0% O.L. 

- VIP 2.5% GAP, 3.1% O.L. 

- VIP 1.5% GAP, 1.0% O.L. 



FIG. 5.13 COMPARISON OF MEASURED AND PREDICTED 
PRESSURE DISTRIBUTION FOR FOSTER'S 
AIRFOIL FLAP COMBINATION 
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Table 2 - RAE 2815 


Flap Configuration 

C L 

meas . 

%»1C. 

C D 

meas. 

C 

D Calc. 

10° flap 
• 0 2 5 C Gap 

1.897 

1.885 

.0222 

.025 

10° flap 
.015C Gap 

1.84 

1.835 

.0203 

.024 

30° flap 
. 020C Gap 

3.40 

3.56 

% .05 - .06 

% . 065 


The infinite swept wing capability of the program is considered in Figures 
5.15 through 5.19. Calculations have been made for comparison with Cumpsty and 
Heads measurements (Ref. 34) on a 61.1° swept wing. Unfortunately the con- 
figuration tested had a large region of separated flow on both the upper and 
lower surfaces. As shown in Figure 5.15 the measured and calculated pressure 
distributions are in reasonable agreement in the forward section of the airfoil 
although the separated flow greatly modifies the pressures in the trailing edge 
region. The predicted streamwise momentum thickness variation is in good agree- 
ment with experiment away from the separation region as shown in Figure 5.16. 

The comparison between predicted and measured values of the angle 3 is good 
only in the region far removed from separation, while the predicted shape 
factor H is in poor agreement with experiment (Figure 5.17). The behavior 
of H is a result of much larger calculated pressure gradients than exist in 
the experimental case. 

Figures 5.18 and 5.19 represent the results of calculations for the RAE 
2815 airfoil flap configuration swept 25 degrees. The pressure distribution 
resulting from 5 iterations of program VI? is shown in Figure 5.18. Also 
included are the calculated aerodynamic lift drag and moment coefficients for 
both the 25 degree and zero degree cases. In this comparison the lift and moment 
coefficients are reduced slightly, while the drag as a result of the increased 
streamwise distance is considerably greater for the swept wing. Flap trailing 
edge streamwise and cross flow velocity profiles are shown in Figure 5,19. For 
this particular case the cross-flow profiles on the flap upper surface are very 
small, a result both of the moderate sweep angle, and the moderate loading on the 
flap. 
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CALCULATIONS VIP 

(RAE 2815 10 FLAP DEFLECTION) 

2.5% GAP, 3.1% O.L. 
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FIG. 5-19 PREDICTED STREAMWISE AND CROSS-FLOW 
VELOCITY PROFILES AT FLAP TRAILING 
EDGE FOR FOSTER'S AIRFOIL FLAP 
CONFIGURATION SWEPT 25 DEGREES. 



PROGRAM LIMITATIONS 


Although it is believed that the calculation procedure and computer program 
is capable of analyzing a wide variety of airfoil configurations, limitations 
both in theoretical methods and due to program structure restrict the range 
of application. These limitations include restriction to: 

Infinite swept wings consisting of at most four elements two of 
which can be slotted flaps. 

Incompressible flow; although the pressure distributions are 
corrected for Mach number effects using Gotherts rule. 

Small regions of separation. Although the source method lends 
itself readily to the development of a separated flow model, the 
current model does not have this capability (as apparent from 
the results of Figure 5.15). If separation is predicted the exist- 
ing approach is to simply extrapolate the source strength to the 
trailing edge of the airfoil. 

Oscillations in lift occurred in early calculations having predicted 
regions of separation. It was finally determined that this was due to un- 
acceptably large source strengths at the trailing edge of the airfoil. 

A numerical experiment demonstrated that monotomically convergent solutions 
can be obtained if the maximum value of the source strength at the trailing 
edge of single element airfoils be limited using the criterion 


where 


v 


ax 


0.115 - 0.01 (ITR-1) 


1 < ITR < ITRMAX 


For configurations having slotted flaps the maximum value of q currently 
allowed in the program is 0.15. In all cases analyzed to date the behavior 
of the lift coefficient has been a monotonic decrease with increasing number 
of iterations. 
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CONCLUSIONS AND RECOMMENDATIONS 


The viscous/potential flow interaction program described in this report 
employs potential flow and boundary layer procedures which are currently 
unique to this method and it is believed that because of this the analysis 
represents a considerable advance on other methods of its type. The program 
in its present form is applicable to a wide variety of problems, in particular, 
to the infinite swept wing case. There is however considerable scope for 
extension to a general three-dimensional wing calculation procedure. 

Specific conclusions regarding each of the major components of the program 
are given in the following paragraphs. 

- The vortex lattice potential flow method developed for this program 
is an accurate numerical approach, adapted for two-dimensions, from a 
general three dimensional lifting potential flow method developed by one of 
the authors. The extension to the three dimensional case is therefore relatively 
straightforward . 


- The integral boundary layer method used for single element airfoils 
and currently on all but the flap upper surfaces in the multi-element mode, 
is quite accurate, as witnessed by the good drag prediction capability. At 
the same time the method uses only a fraction of a second of computer time 
per boundary layer development. 

The inclusion of curvature and normal pressure gradient effects in the 
finite difference boundary layer method enables complicated boundary layer 
flows to be represented more accurately than is possible by other procedures 
now available. It is believed that the improved physical representation of 
the flow over slotted flaps is responsible for the greater degree of agreement 
with experiment than is achieved by other methods. 

— The use of sources to represent the displacement effects of the boundary 
layer on the potential flow, while not necessarily more accurate a procedure 
than the direct employment of the displacement thickness, has two distinct 
advantages. The first advantage is in the computational superiority of such a 
procedure. The influence coefficient matrix need by inverted only once, with 
succeeding iterations requiring only matrix multiplication. If one is ever 
to contemplate a viscous/potential flow interaction program applied to general 
three dimensional airplane configurations such a procedure will be almost 
mandatory. The second advantage is related to the modelling of separated flow 
regions in a potential flow analysis by distributed sources. It has been 
demonstrated in the literature that such an approach is possible, and it is 
recommended that one of the first extensions of the present program should 
be to include separation effects. The program would then be capable of 
predicting C as a function of Reynolds number, 
max 
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Another useful extension to the program would be the inclusion of 
compressibility effects in the boundary layer development. Even at low free 
stream Mach numbers the high suction peaks experienced by a slot or main element 
of a high lift system can lead to compressibility problems. 


With some modification the finite difference boundary layer program 
module can be used to predict the effect of tangential blowing or boundary 
layer suction in the overall context of a viscous/potential flow interaction 
method for high lift systems. 
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APPENDIX I 


POTENTIAL FLOW THEORY 


The potential flow theory is used to derive the influence of constant 
and linear distributions of sources and vorticity on planar two dimensional 
surfaces. Consider an elementary line source or line vortex located at a point 
£ on the x axis and perpendicular to the x, z plane. In incompressible flow, 
the magnitude of the velocity induced by either singularity at an arbitrary 
point P(x,z) is given by: 


where 


V = 


1 

2?rd 


( 1 ) 


d = [(x-5) 2 + z 2 ] 


L /2 


The geometry is illustrated by the following sketch: 



For an elementary source, the velocity V is in the direction of the line 
joining £ and P, while for an elementary vortex, the velocity V is 
perpendicular to this line. The horizontal and vertical components of velocity 
corresponding to the line source or vortex are given by: 


w = u “ Vsin0 
s v 



( 2 ) 
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u = - W = VCOS0 - (3) 

s V 2,d 2 

The contribution of a constant distribution of sources or vortices along the 
x axis is obtained by integrating equations (2) and (3) from 0 to c, 


w 

s 


u 

s 


= u 


-£/ 


dg 

2 2 
(x~U + z Z 


1 [ -1 Z 

= -r— tan 

2tt x~c 


tan 


-1 


= -w 


2ir J 


(x-g)dg 

2 2 
(x-5) + z Z 




/(x -c) 2 +z 2 

Vx 2 + z 2 


(A) 


(5) 


The effec ts of compressibility may be obtained by applying Gothert’s Rule, 
with B = /1 -m2~ 


and 


w 


s 


u 

V 


2tt 



Bz 

x-c 


tan 


-1 



( 6 ) 


u 

s 


i /(x-c ) 2 + e 2 z 2 

2 "8 08 V * 2 + e 2 z 2 


(7) 


ft 2 

w = - P u 
v s 


_6_ 

2ir 


log 


/(x-c) 2 + e 2 z 2 
V x 2 + e 2 z 2 


( 8 ) 


For a source or vortex distribution varying linearly with x, with zero strength 
at the origin and unit strength at x = c. 


w 

s 


= u 

V 


z 

2ttc 


/ 


m 

2 2 
(x-C) + z 
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"2tt 



u 

s 


-w 

V 


1 

2ttc 


f S(x-C) d£ 

J r j. 2 

o (x-£) + z 


^1 

2tt 




£ 

x 


-1 

tan 


z 

x-c 


+ C 108 



( 10 ) 


Compressibility effects are obtained as before, i.e. , by multiplying z and 
w v by and dividing u s by 8. 

A source or vortex distribution having unit strength at the origin, and 
zero strength at x = c, can be obtained by subtracting the previously derived 
linearly varying distributions from the constant distributions. 
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APPENDIX II 


SOLUTION OF BOUNDARY CONDITION EQUATIONS 


For single or multi-element airfoils, the boundary condition equations can 
be solved by direct inversion. For multi-element airfoils, use can also be 
made of a rapidly convergent iteration scheme reported in Reference (35) . In 
this method the matrix is subdivided into smaller partitions, or blocks, with 
each block representing the influence of one element of the multi-element air- 
foil. The diagonal blocks represent the influence of the elements on themselves, 
the off-diagonal blocks represent the interference of one element on the others. 
The order of any block is restricted to 60, the maximum number of panels on the 
upper and lower surface of the element. 

The initial iteration calculates the source and vortex strengths 
corresponding to each block in isolation. For this step, only the diagonal 
blocks are present in the aerodynamic matrix. Once the initial approximation 
to the source and vortex strengths is determined, the interference effect of each 
block on all the others is calculated by matrix multiplication. The incremental 
normal velocities obtained are subtracted from the normal velocities specified 
by the boundary conditions. This process is repeated 15 to 20 times, or until 
the residual interference velocities are small enough to ensure that convergence 
has occurred. 

The procedure is illustrated below for an aerodynamic matrix consisting 
of nine blocks. The unknown singularity strengths are designated the 
specified normal velocities C^. ^ 


To solve 


where 



A 11 

A 12 

A 13 

A 

A 

A_ 

21 

22 

23 

A 

A „ 

A„„ 

31 

32 

33 
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Put 


A = D + E 


A 11 ° 


22 


0 

0 

A 


33 


Therefore [D + E] {y} = {C} 

or {y} = D -1 [C - E {y}] 

First approximation: 

{y } -1 = D _1 {C} 


*12 


21 


A 31 A 32 


Calculate 


AC 1 = E {y } 1 


Second approximation: 


{y } 2 = D 1 {C - AC 1 } 


til 

Similarly, k approximation: 


Note that 


{y } k = D _1 {C - AC k 1 } 


-1 


aT 1 o o 

0 A 22 0 

0 0 


si 
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APPENDIX III 


PROGRAM MACRO FLOW CHARTS 


The general program overlay structure of each overlay is described 
briefly in the following paragraphs and flow charts. * 

OVERLAY (0,0) Program VIP 

Program VIP controls the entire computer program (Figure Al) . All 
primary overlays are called from VIP. Overlay (0,0) also contains other 
subroutines which are commonly used in the other overlays. 

OVERLAY (1,0) Program POTFLOW 

Program POTFLOW controls the lofting of the configuration, the calculation 
of flap surface curvature, the calculation of the potential flow pressures, as 
well as the calculation of the lift and moment coefficients. 


OVERLAY (2,0) Program IBL 

Program IBL controls the integral boundary layer analysis from the 
calculation of initial conditions along a stagnation line to the turbulent 
boundary layer analysis. The program logic flow is shown in Figure A3. 

OVERLAY (3,0) Program INSPAN 

Program INSPAN controls the infinite swept wing finite difference boundary 
layer analysis. The overlay calls two secondary overlays, OVERLAY (3 >1) Program 
BOUNDRY (Figure A4) and OVERLAY (3 ,2) , Program DEVELOP (Figure A5) , Program 
BOUNDRY initializes the grid network normal to the flap surface upon which the 
finite difference method is applied. Normal chord and spanwise velocity profiles 
are initialized in preparation for the analysis. 

Program DEVELOP controls the actual calculation procedure used in determining 
the downstream development of the boundary layer. 


OVERLAY (4,0) Program FELDPT 

Program FELDPT calculates the off-body pressure distributions P(x,y) for 
input to INSPAN. 

* Program input/output discription and a complete program listing is given in 
Reference 36. 
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FIG. A1 OVERLAY (0,0) PROGRAM VIP 





FIG. A -2 OVERLAY (1,0) PROGRAM POTFLOW 
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FIG. A3 OVERLAY (2,0) PROGRAM IBL 


84 










FIG. A - 4 OVERLAY (3, 1) PROGRAM BOUNDARY 
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FIG. A - 5 OVERLAY (3, 2) PROGRAM DEVELOP 


NASA - Langley, 1974 
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